Stability of three-sublattice order in S = 1 bilinear-biquadratic Heisenberg Model on 

anisotropic triangular lattices 
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The S — 1 bilinear-biquadratic Heisenberg model on anisotropic triangular lattices is investigated 
by several complementary methods. Our focus is on the stability of the three-sublattice spin nematic 
state against spatial anisotropy. We find that, deviated from the case of isotropic triangular lattice, 
quantum fluctuations enhance and the three-sublattice spin nematic order is reduced. In the limit 
of weakly coupling chains, by mapping the systems to an effective one-dimensional model, we show 
that the three-sublattice spin nematic order develops at infinitesimal interchain coupling. Our 
results provide a complete picture for smooth crossover from the triangular-lattice case to both the 
square-lattice and the one-dimensional limits. 

PACS numbers: 75.10.Jm, 75.10.Kt 



I. INTRODUCTION 

Spin nematic states are the states of quantum spin 
systems in which no spin-dipolar ordering exists, but 
spin-rotation symmetry is spontaneously broken due to 
the appearance of spin-quadrupolar order. [1] Promi- 
nent examples for the existence of such spin nematic 
phases include the spin-1 bilinear-biquadratic (BLBQ) 
model. [1, 2] Interest in quantum states with spin nematic 
order has been raised recently by experimental findings 
in NiGa2S4, which is an insulating quantum magnet with 
spin-1 Ni 2+ ions living on a triangular lattice. [3] This 
system is found to be in a gaplcss ground state with- 
out spin-dipolar ordering. It has been suggested that 
this compound can be considered as a physical real- 
ization of the BLBQ model on a triangular lattice and 
the candidate ground state is characterized by a three- 
sublattice spin nematic order. [4, 5] The observed gap- 
less excitation spectrum thus corresponds to the Nambu- 
Goldstone modes associated with spontaneous breaking 
of spin-rotation symmetry. 

While consensus has been reached for the ground states 
of the BLBQ model on a triangular lattice, [1] physics for 
spatially anisotropic models has not yet been addressed. 
Here we consider the spin-1 BLBQ model on anisotropic 
triangular lattices [see Fig. 1(a)] defined by the Hamilto- 
nian, 

H = Ji [ cosd s i ■ s j + sin0 ( s i ' s i) 2 ] 

(id) 

+J 2 [ cos0 S * ' S i + sin6> ( Si ' S J') 2 ] ' w 

where Si's are spin-1 operators. We use the notations 
and ((«', j)) to denote the nearest-neighbour bonds 
and the bonds along only one of the diagonals, respec- 
tively. Ji and Ji are the coupling strengths on the corre- 
sponding bonds. The relative strength of the linear and 
the biquadratic couplings is parameterized by 9. Here 
a = J2/J1 defines the extent of spatial anisotropy. As 
the anisotropy a increases from zero, the model changes 
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FIG. 1: (Color online) (a) Illustration of the anisotropic tri- 
angular lattice in a square topology and the schematic repre- 
sentation of three-sublattice spin nematic order. Two groups 
of interactions, Ji and J2, are denoted by solid and dashed 
links, respectively. Here the three mutually orthogonal vec- 
tors di in the mean-field analysis (see Sec. II A) are associated 
with three different colors, (b) Brillouin zone of square lat- 
tice. The reduced Brillouin zone for the three-sublattice order 
is enclosed by dashed lines. The k path used in Fig. 4 is de- 
fined as follows: Y : (0,0), M : (ir/3,ir/3), K : (-7r/9, 7n/9), 
and Q : (-5tt/9, 5tt/9). 



from the square lattice to the isotropic triangular lattice, 
and eventually to decoupled chains. We concentrate on 
the parameter region of Ji, J 2 > and 7r/4 < 9 < it/2, 
where the ground states with three-sublattice spin ne- 
matic order are expected. 

The model in Eq. (1) includes several limiting cases, in 
which the ground states are known: 

(i) At Ji = 0, the anisotropic model becomes a set 
of decoupled one-dimensional (ID) BLBQ spin chains, in 
which each spin interacts with two neighbors only (i.e., 
the coordination number z = 2). For each spin chain 
with 7r/4 < 9 < 7r/2, the system is found to be in an ex- 
tended critical phase with soft modes at momenta k = 0, 
±2tt/3. [6] Away from the SU(3) point (9 = tt/4), this 
phase develops dominant antiferro-quadrupolar correla- 
tions with a period of three lattice units (i.e., almost 
"trimerized" ground state). [7, 8] 
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(ii) At J i = J2, the model in Eq. (1) is equivalent to 
an isotropic triangular-lattice model with z = 6. The 
ground state for 7r/4<#<7r/2is shown to poss a three- 
sublattice spin ncmatic order, [4, 5] where the nematic 
directors on the three sublattices A, B, and C of the 
triangular lattice are orthogonal to each other (say, along 
x, y, and z, respectively). The schematic representation 
of this order is shown in Fig. 1(a). 

(hi) At J2 = 0, our model reduces to a square-lattice 
model with z = 4. It is established only recently that the 
ground state for tt/4 < 9 < tt/2 develops an unexpected 
three-sublattice spin ncmatic order as a consequence of 
a subtle quantum order- by-disorder mechanism. [9-11] 

In this paper, the spatially anisotropic BLBQ model 
in Eq. (1) is investigated. We pay our attention to the 
effect of spatial anisotropy on the stability of the three- 
sublattice spin ncmatic state in this model. Our main 
results for generic cases of anisotropy are based on the 
linear flavor-wave (LFW) theory, [1, 2, 12, 13] which 
has been applied to the triangular-lattice as well as the 
square-lattice cases with success. [4, 5, 9-11] We find 
that three-sublattice spin nematic order is most robust 
in the case of isotropic triangular lattice with anisotropy 
a = J2/J1 = 1- As deviated from this a — 1 case, quan- 
tum fluctuations enhance and the order is reduced. This 
behavior is reasonable, since the coordination number z 
is decreased both in the a — > (square-lattice limit) and 
the a — > 00 (decoupled-chain limit) cases, and stronger 
quantum fluctuations are thus allowed. In order to ad- 
dress the validity of the LFW theory, we have performed 
exact diagonalizations (ED) on lattices of small sizes. By 
comparing our LFW predictions specialized to finite-size 
systems with the numerical results, we find that quantum 
fluctuations obtained by the LFW theory are overesti- 
mated, especially in both a — > and a — > 00 limits. Thus 
the stability region of the three-sublattice state could be 
larger than the LFW predictions. Since previous numeri- 
cal investigations, [9-11] have shown nonzero order in the 
square-lattice case at 9 = 7r/4, one may expect that the 
three-sublattice order could persist down to the a = 
limit in the whole region of7r/4 < 9 < tt/2. In the oppo- 
site large-a limit, the status is much less clear. Because 
there is no true long-range order at the decoupled-chain 
limit (J1/J2 = 1/cx = 0), [6-8] an interesting issue is 
whether a nonzero interchain coupling is necessary or 
not for the appearance of two-dimensional (2D) three- 
sublattice order. By mapping from the system of weakly 
coupled chains (J1/J2 -C 1) to an effective ID model, 
we show that the critical value of interchain coupling is 
{Ji/J2)c = for all 7t/4 < 9 < tt/2. In other words, the 
transition from the 2D three-sublattice phase to the ID 
"trimcrized" critical phase [6-8] should occur at infinite 
a. 

The rest part of the paper is organized as follows. 
Generic cases of anisotropy are discussed in Sec. II, where 
details of the LFW analysis are presented in Sec. II A and 
the comparison between the finite-size LFW and the ED 
results is made in Sec. II B. The case in the decoupled- 



chain limit is explored in Sec. Ill through field-theoretical 
approach as well as ED calculations. The last section is 
devoted to our conclusions. 



II. GENERIC ANISOTROPY 

A. linear flavor-wave analysis 

The LFW theory starts from representing the model in 
Eq. (1) in terms of three-flavor Schwinger bosons a^ a un- 
der the local constraint ^2 a a\ a ai >a = 1. [1, 2, 4, 5, 9-13] 
The Schwinger bosons a\ a (with a = x, y, z) create three 
time-reversal-invariant local basis states, \x) = ~^{\ s z = 
1) - \s z = -1)), \y) = -±(\s z = 1) + \s z = -1)), and 
\z) = —i \s z =0). In terms of these bosons, the spin 
operators become Sf = —i Y^r e a| g 7 aj aCii, 7 . We denote 
di as the local ordering vector and let {dj, e^, f^} forming 
a local orthonormal basis. A generic local quantum state 
can be represented by the linear combination of the three 
basis states {|dj) , |ej) , |fi)}. Let a\, b\ and c] represent- 
ing the Schwinger boson operators which create the local 
states {|dj) , |e^) , \fi)} out of the Schwinger boson vac- 
uum. They are related to the operators Oi ;Q , through the 
relation a; Q = + ej, Q bi + fi^ a Ci. Within the LFW 

analysis, we solve the local constraint a\di + b\bi + c\ci = 

1 by replacing a; = y 1 



1 and thus 



a-i.a ~ di,a + P-i.a bi + fi. a Ci. It has be shown that, for 
7r/4 < 9 < n/2, the mean-field energy of the nearest- 
neighbor bond is minimized when the d^ vectors are mu- 
tually orthogonal. [1] In the case of isotropic triangular 
lattice, these d^ vectors are given by the unit vectors 
along the x, y, and z directions on the three sublattices 
[see Fig. 1(a)]. Employing this mean-field condition and 
the approximate expression for a at i, the model in Eq. (1) 
reduces to the following quadratic LFW Hamiltonian (up 
to a constant term), 

#lfw = 2 J2 [eo^bk + 4°k) + (Ak^kc-k + hue.) 



(2) 



Here the values of k run over the first Brillouin zone of 
the square lattice and 



£o = (^1 + y) sin( 



A k = 



cos 9 



,h{e lk * + e lk y) + J 2 e - t{k * +k «A , (3) 



</> k = (tan6>- l)A k . 

For k ^ 0,±k with k = (2tt/3, 2tt/3), the resulting 
quadratic bosonic Hamiltonian can be diagonalized by 
the Boguliubov transformation, and the corresponding 
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excitation spectrums are given by 



wi, a (k) =2 



4 + 



±j2\^(2el-\A k \ 



1/2 



(4) 



The modes for k = 0, ±k cannot be diagonalized in this 
way, because the Boguliubov transformation becomes 
singular here. As discussed in the finite-size spin-wave 
theory for spin-1/2 Heisenberg model, [14, 15] these sin- 
gular modes have no contribution to the ground-state 
energy, while removal of these modes is required in the 
computation of order parameter. 

Some general features of the LFW excitation spec- 
trums are described below. At the SU(3) point of 9 = 
7r/4, we have </>k = and thus these two excitation modes 
become degenerate in energy. Away from this special 
point, wi(k) gives a gapped mode, while w 2 (k) is gaplcss 
and has nodes at k = 0, ±ko. We remind that the primi- 
tive unit cell of three-sublattice states contains three lat- 
tice sites as its basis and therefore its size becomes three 
times larger. As a consequence, the original Brillouin 
zone of square lattice behaves as an extended Brillouin 
zone [as shown in Fig. 1 (b)] , such that each branch of ex- 
citations in Eq. (4) becomes three-fold degenerate within 
the original Brillouin zone. As a simple check for our 
derivations, we point out that the obtained dispersions 
at J 2 = do reduce to those in Ref. 9 for the square- 
lattice case. For the case of isotropic triangular lattice 
(Ji = J 2 ), they are equivalent to the results in Ref. 4. 

To determine the stability region of the three- 
sublattice states, a suitable order parameter should be 
measured. In the spin nematic state, spin rotational 
symmetry is spontaneously broken, though time rever- 
sal symmetry is preserved. In such a state, average 
magnetic moment must vanish ((S) = 0). Nevertheless, 
quadrupole order can appear, which is characterized by a 
nonzero expectation value of the symmetric and traceless 
rank-2 tensor operator 



Q 



sf + S?S? 



(5) 



Here Sf is the a component of spin-1 operator at site i 
and 8 a P is the Kroneker delta symbol. In terms of the 
local ordering vector d^, the expectation value of this 
quadrupole operator can be written as 



f) = -q( d?df - U afS 



(6) 



where the constant value of q describes the magnitude of 
the quadrupolar ordering. For both cases of the trian- 
gular and the square lattices, [4, 9] the vectors of the 
three-sublattice states point along three orthogonal direc- 
tions in three different sublatticcs, as shown in Fig. 1(a). 
From Eqs. (5) and (6), we have (from now on, summation 
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FIG. 2: (Color online) Effect of anisotropy a = J2/J1 on 
the quadrupole order parameter q within the LFW theory for 
various values of 9. The inset shows the details around the 
region of a = 0. 



is implied over the repeated Greek indices) 

q = -l{Qf)d<*4 = l-|((S i -d,) 2 ). (7) 

At classical level, q = 1 because (S, • dj)|dj) = 0. 
Within the LFW theory, q can be expressed by 



q = 1 - - (An„) 



(8) 



with 



(An a > = 1 - 1 £ (aja.) = 1 £ (b\ bi + cfo) (9) 



being the deviation of the number density for the 
Schwinger boson from its classical value of one. Here 
N is the total number of lattice sites. We note that this 
expression of q is nothing but the local moment defined 
in Eq. (24) of Ref. 11. From this expression, it is obvious 
that the quantum correction for q comes from the non- 
vanishing contribution of (An a ). When (An Q ) increases 
to 2/3, q vanishes. This gives a phase transition out of 
the three-sublattice states within the LFW theory. The 
explicit expression of (An a ) is given by 



(An n 



- E 

N ^ 



£o + \<h\ , £ o - 1 0k I 



N ^ V wi(k) w 2 (k) 



1 



(10) 

Here the singular modes at k = 0, ±ko are excluded from 
this summation. Note that, for the case of square lattice 
(J 2 = 0) and at the SU(3) point (6 = tt/4), Eq. (10) 

reduces to (An a ) = i £k#o,±k„ ( ^q^p ~ X ) with 

7k = cos[(k x — k y )/2], and reproduces the previous re- 
sult (see Eq. (22) of Ref. 11). The general behavior of 
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FIG. 3: (Color online) Phase diagram for the spatially 
anisotropic 5 = 1 BLBQ model in Eq. (1) determined by 
the LFW theory. The insets show the details of two phase 
boundaries around 8 = 7r/4. 
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the order parameter q as functions of a for distinct val- 
ues of 9 is shown in Fig. 2. It is seen clearly that the 
three-sublattice nematic order is most robust at a = 1. 
Far away from this point of isotropic triangular lattice, 
quantum fluctuations arising from the flavor- wave excita- 
tions become more stronger, and they destroy eventually 
the quadrupolar ordering in both limits of a — > and 
a. — y oo. Exploiting Eqs. (8) and (10) and employing 
the condition q = as the criterion for the transitions 
out of the three-sublattice states, we can establish the 
phase boundaries of the three-sublattice states as shown 
in Fig. 3. We find that, in general, the lower transition 
points are nonzero and the upper ones are large but finite. 
The stability region of the three-sublattice order is largely 
reduced as 9 approaches the SU(3) point (9 = tt/4). It 
implies that there exist more low-lying excitations and 
thus larger quantum fluctuations as 9 gets closer to 7r/4. 

As seen from the expression of Eq. (10), the gapless 
mode W2(k) should make a dominant contribution in re- 
ducing the three-sublattice order. To have a better un- 
derstanding of the enhancement of quantum fluctuations 
both in the square-lattice limit (a —> 0) and the quasi-lD 
limit (a 3> 1), it is instructive to examine the softening 
behavior of this excitation mode more closely. The flavor- 
wave dispersions of the gapless branch W2(k) along the 
path defined in Fig. 1(b) for various values of anisotropy 
a = J2/J1 are shown in Fig. 4. It can be seen that, as 
system approaches the square-lattice limit (a — > 0), the 
flavor-wave velocity at the T point (defined by the slope 
of the dispersion relation) decreases to zero. Therefore, 
the excitation modes along the T-M line (i.e., line of 
k x = k y ) become zero-energy modes eventually. On the 
other hand, in the limit of decoupled chains {a — > 00), 
the excitation energies in the whole Brillouin zone are 
softened. Within the LFW theory, the disappearance of 
three-sublattice order in both limits of a — ¥ and a — > 00 
can be explained by such softening in energy. 



FIG. 4: (Color online) Dispersion relation of gapless branch 
cj2(k) of the flavor- wave excitation for different anisotropy 
parameters a = J2/J1 at 9 — 0.3ir. Note that energies are 
measured in unit of Ji for a < 1 (upper panel), while they 
are measured in unit of J2 for a > 1 (lower panel). 



Near the node, say, at k = 0, analytic expressions can 
be derived. By expanding Ak in Eq. (3) near k = 0, 
we can show that the gapless flavor-wave mode behaves 

as W2(k)/Ji w c 2 + k\ + c^k 2 . with k± = ^{k x ± k y ), 

c + = ^^/sin(26»)a, and c_ = ^^sin(26»)(a + 2). For 
a — > 0, we have c + — > while c_ remaining finite. The 
outcome of c + = for a = gives a nodal line in the 
excitation spectrum along the k x = k y direction (i.e., T- 
AI line), as observed in Fig. 4. Within the LFW theory, 
these soft modes play a significant role in the destruction 
of the nematic order in the square-lattice limit, as no- 
ticed in the previous investigations. [9, 11] On the other 
hand, in the extreme anisotropic quasi-lD limit ( J\ — > 
or a — >• 00 ) , we have c_ /c + = 1/3 for all values of 9. That 
is, this ratio of the flavor-wave velocities does not goes to 
zero in the quasi-lD limit. Instead, the whole spectrum, 
in unit of the diagonal coupling J2, becomes nearly flat 
in the entire Brillouin zone, as seen from Fig. 4. There- 
fore, the associated quantum fluctuations become more 
and more significant and finally the 2D nematic order 
ceases to exist for a being large enough. We stress that 
the mode softening of the flavor-wave excitations in the 
quasi- ID limit is quite different from what one got for 
the spin- wave excitations in spatially anisotropic spin-1/2 
antifcrromagnetic Heisenberg models on cither triangu- 
lar [16] or square lattices. [17, 18] In these spin-1/2 cases, 
one can show that, as the interchain couplings approach 
zero, only the spin-wave velocity c± for the excitation 
transverse to the chains will vanish, but the spin-wave 
velocity cm for the excitation along the chains will remain 
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finite. Thus the whole spectrum, in unit of the intrachain 
coupling, never becomes nearly flat in the whole Brillouin 
zone, and the ratio of these two spin- wave velocities ci/cy 
does go to zero in the quasi-lD limit. 

We remind that the LFW analysis is valid only when 
quantum fluctuations arc weak (i.e., only when (b\bi) + 
(c\ci) <C 1) because the local constraint, aja 2 ; + b\bi + 
c\ci = 1, is considered only approximatively. Therefore, 
we should be cautious with the LFW results about the 
phase boundaries shown in Fig. 3, since large quantum 
fluctuations are expected near the transition points. To 
examine the validity of the LFW predictions, comparison 
with exact results is necessary. 
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B. exact diagonalizetion 

In this subsection, we perform ED calculations for the 
Hamiltonian in Eq. (1) on small clusters and compare the 
results with those obtained by the LFW analysis on the 
same clusters. 

We remind that there exist subtleties in making care- 
ful comparison of order parameter between symmetry- 
breaking solutions (say, LFW results) and symmetry- 
nonbreaking ones (say, ED findings). Such an observation 
has been put forward in Ref. 19 in concern with magnetic 
ordering in spin- 1/2 Hcisenberg antiferromagnets on a 
triangular lattice. To uncover the long-range order on 
lattices of small sizes, a proper quantity has to be mea- 
sured in ED calculations. Here the squared quadrupolc 
moment Q 2 in a given sublattice (say, A sublattice) is 
considered, 




As mentioned before, the Einstein summation conven- 
tion for the repeated Greek indices is assumed. The sig- 
nature of three-sublattice order will be manifested as a 
macroscopic value of Q. To obtain a order parameter 
that is normalized to 1 in the absence of quantum fluc- 
tuations, the sublattice quadrupolc moment Q should be 
divided by a size-dependent normalization factor. For de- 
riving the correct normalization factor, it should be kept 
in mind that the sublattice quadrupolc moment cannot 
be treated as a classical quantity. For example, it can 
be shown that there exists an exact operator identity: 
Qf 3 Q"' 3 = 5/3 on a given site i. On the other hand, 
when i and j denote different sites of the same sublat- 
tice, Q°/®\ = 2/3 for fully aligned classical ordered 

state. From these observations, the maximum quantum 
value of Q 2 can be shown to be (2N 2 /27)(1 + 9/2N) for 
systems of N sites. Thus a valid definition of the order 
parameter would be 

j 27Q 2 

9= V^oW (12) 



FIG. 5: (Color online) Comparison of order parameter q be- 
tween the ED and the LFW results for different lattice sites as 
anisotropy a ia varied. Here, = arctan(2). The ED data are 
obtained on lattices with sites iV = 9 (triangles) and N = 12 
(squares). The LFW results of the corresponding sizes are 
shown by solid and dashed lines, respectively. The lowest 
curve refers to the infinite-size LFW results (dotted line). 



Now q saturates at one in the classical state and should 
be decreased by quantum fluctuations in the quantum 
ground state. We stress that, for careful comparison be- 
tween the ED and the LFW results for small values of N, 
it is important to use the correct normalization factor in 
the definition of the order parameter q. 

The comparison between the ED and the LFW results 
is shown in Fig. 5 for 8 = arctan(2). The ED data are 
calculated by using Eqs. (11) and (12). On the other 
hand, the LFW results for systems of finite sizes arc 
obtained from Eqs. (8) and (10) by summing the mo- 
menta (except k = 0, ±kn) determined by the clusters 
in ED calculations. We find that, around the isotropic 
point (a = 1), two sets of results do not differ by large 
amounts. Good agreement can persist even down to the 
a = limit for the special case of A" = 9. This indi- 
cates that the LFW theory does in general provide good 
approximation around a = 1. However, serious reduc- 
tion in the LFW results of q as compared to the ED 
ones is observed both when a -C 1 and a>l, This im- 
plies that quantum fluctuations are significantly overesti- 
mated in the LFW analysis in both of the square-lattice 
and the quasi-lD limits. In other words, the softening 
of the flavor-wave excitations in both limits (see Fig. 4) 
should be exaggerated. Thus one has to go beyond the 
LFW approximation to find improvements on the exci- 
tation spectrums. According to previous numerical in- 
vestigations, [9-11] where nonzero order was reported 
in the square-lattice case at 9 = it/ 4, the lower phase 
boundary obtained within the LFW theory (see Fig. 3) 
may be illusive. Instead, the three-sublattice order may 
persist down to the a = limit in the whole region of 
tt/4 < 6 < it/2. 

The above comparison suggests as well that the upper 
phase boundary may take much larger values than the 
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ones estimated by the LFW theory. To achieve the true 
values of the transition points in the large a limit, it is 
instructive to analyze the model from its quasi-lD limit. 
Because there is no true long-range order in strictly ID 
models, our main concern is to show whether a nonzero 
interchain coupling is necessary or not to establish the 
2D order. This is what we shall do in the next section. 



III. WEAKLY-COUPLED-CHAIN LIMIT 

When Ji <C J2, the system described by Eq. (1) re- 
duces to weakly coupled ID chains with intrachain cou- 
pling strength J2 and weak interchain coupling strength 
J\. Taking advantage of conventional mean- field treat- 
ment for the interchain coupling, [17, 18, 20, 21] our 
quasi-lD systems can be transformed into effective single- 
chain problems. 

The desired effective single-chain model can be de- 
rived in the following way. Using the definition of the 
quadrupole operator in Eq. (5), the part of Hamiltonian 
with the interchain coupling J\ can be rewritten as 

H x = Jocose- ^) £ Si • Sj + Ji sin0£ Qf Qf 



(13) 

by dropping some constant terms. Again, the repeated 
Greek indices imply the Einstein summation convention. 
Assuming the three-sublattice quadrupolar ordering, in 
which (Sj) = but (Qf) is nonzero, Hi can be approx- 
imated by an on-site Hamiltonian, 



#1 



Jisin0£(Q 



a/3 i 
n+m/ 



(14) 



where m runs over all neighbors on the J\ bonds for 
the n-th site along a given single chain. Substituting 
suitable mean-field expression for (Q"+ m ), an effective 
ID Hamiltonian for the model in Eq. (1) can be written 



J2 £ [cos 9 S„ • S„ + i + sin 9 (S 



-H x . 

(15) 

This describes a ID BLBQ chain in a self-consistent ex- 
ternal field triggering three-sublattice quadrupolar order- 
ing. For convenience, we set J2 = 1 as the energy unit 
in this section. We note that the self-consistent field is 
proportional to Jisin#. It is thus expected that, for a 
given anisotropy (i.e., for a fixed value of Ji), the result- 
ing order will be stronger as 9 gets closer to tt/2. This 
observation is consistent with our LFW results, where 
the stability region of the three-sublattice state is pushed 
toward larger a = J2/J1 as 9 — > tt/2. 

In the following, both analytical and numerical tech- 
niques are exploited to determine the critical interchain 
coupling Ji. c for the emergence of 2D three-sublattice 
order. 



A. scaling analysis 

Employing the mean- field expression of (Qf) in 
Eq. (6), the on-site part of Eq. (14) for the effective ID 
Hamiltonian becomes 



i-?(s, 



(16) 



where the self-consistent field conjugate to the operator 
for the order parameter q in Eq. (7) is defined by h = 
(4/3)gJisin0. 

For Ji <C 1, the effective ID Hamiltonian in Eq. (15) 
can be considered as a ID BLBQ model in a weak exter- 
nal field h. Thus it should be valid to treat the effect of 
h as a perturbation. It is known from Ref. 7 that, near 
9 = 7r/4, the ID BLBQ spin chain can be described by 
an SU(3)i Wess-Zumino-Witten conformal field theory 
perturbed by some marginally irrelevant current-current 
interactions. To see whether the three-sublattice order 
can be induced by vanishing self-consistent field h or 
not, we need only to calculate the scaling dimension Ah 
of the corresponding operator in the on-site term and 
then determine its relevancy. If Ah < 2, the on-site 
term provides a relevant perturbation and thus the three- 
sublattice order will be induced by an infinitesimal h. 
Otherwise, the on-site term becomes irrelevant and the 
three-sublattice order can be established only when h ex- 
ceeds a nonzero critical field h c . In this latter case, we 
need to determine h c numerically by solving this model 
explicitly. 

In terms of the field-theoretical variables discussed in 
Ref. 7, the operator in the on-site term of Eq. (16) takes 
the form of the primary fields of an SU(^) Wess-Zumino- 
Wittcn model for v = 3. That is, the scaling dimension 
A/j of our operator is nothing but that of those primary 
fields, which is equal to 1 — 1/^ = 2/3 according to the 
analysis in Ref. 7. Because of A^ < 2, the perturbation 
caused by the self-consistent field is strongly relevant. 
Since the above scaling argument is essentially indepen- 
dent of the value of 9, we claim that h c — and thus 
Ji.c = for 7r/4 < 9 < tt/2, even though the field the- 
ory in Ref. 7 is derived for 9 close to 7r/4. This implies 
that, for original quasi-lD anisotropic BLBQ model, true 
phase transitions out of the three-sublattice states actu- 
ally occur at infinite a, rather than at large but finite a c 
as obtained in the LFW analysis. 

We can go one step further to establish a nonper- 
turbativc relation between the order parameter q and 
the weak interchain couplings J\ by making use of the 
field-theoretic approach. According to the standard scal- 
ing argument, [22] the order parameter q induced by 
the perturbation of self-consistent field h scales as q oc 
h Ah '( 2 ~ Ah K Combined with the self-consistency rela- 
tion, h = (4/3)gJi sin 6, we get q oc (Ji sm9) A '^^ 1 ~ A ^ . 
Since A^ = 2/3, we conclude that q oc Jisin#. Inter- 
estingly, this result coincides with the one that is ex- 
pected naively from the perturbation theory for original 
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Hi without taking mean-field approximation. This seems 
to indicate that it is possible to study this model in its 
quasi- ID limit directly from perturbation theory for Hi. 
Instead of pursuing along this direction, we shall deter- 
mine the phase boundary by the numerical ED method 
below. 



B. exact diagonalization 

In this subsection, the critical values of the interchain 
coupling Ji are estimated by the ED method. Here we 
follow the treatment in Rcf. 17 for quasi- ID Heisenberg 
antiferromagnets. Our ED results provide numerical evi- 
dences in supporting the above conclusions based on scal- 
ing arguments. 

For the sake of ED calculations, we assume here that 
only the zz component of the expectation value (Q"' 3 ) is 
nonzero. Thus the effective ID Hamiltonian in Eq. (15) 
has still spin-rotation symmetry in the z direction and 
the total z-component spin remains a conserved quan- 
tity. This reduces much computational effort and thus 
calculations for systems of large sizes become available. 
In consistent with Eq. (6), the explicit form of (Q"^) is 
taken to be (Qf ) = ((S"|) 2 ) - | = -|gcos(Q • r,-). Sub- 
stituting the present mean-field solution to Eq. (14), the 
on-site part of the effective ID Hamiltonian becomes 

Hi=h^2cos( Y n)(S z n f . (17) 

n 

Here the self-consistent field is again given by h = 
(4/3)gJisin0. 

By taking h as a free parameter, we diagonalize nu- 
merically the effective single-chain model up to system 
length L = 18. For the present single-chain problem, the 
order parameter for spin quadrupole ordering becomes 

g = -|Xy*"<0 . (18) 

n 

This expression is compatible with the form of 2D order 
parameter used in this subsection. The results of q as 
function of h for several O's with L = 18 are presented in 
Fig. 6(a). The susceptibility x = (dq/dh)\ h=0 can then 
be evaluated from the slope of the linear fit as shown in 
the inset of this figure. Within the present mean-field 
approach, to have a nonzero solution of q, the slope x °f 
the tangent line around h = must be larger than that of 
the straight line, q = /i/[(4/3) Ji sin#], given by the self- 
consistent relation. That is, long-range order appears 
only when x > l/[(4/3)Ji sin#]. This requirement leads 
to a critical value Ji tC of the interchain coupling for a 
given length L, 

1 

Jl ' c = (4/3) X sin0 ' (19) 

The size dependence of Ji jC for various O's is shown in 
Fig. 6(b). As seen from this figure, size dependence of 
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FIG. 6: (Color online) (a) Order parameter q as function of 
self-consistent field h for various 9's with L = 18. Lines are 
guide to eyes. Inset: linear fit around h = region, (b) 
Critical value Ji, c of the interchain coupling as function of 
1/L. Lines show the extrapolations in the thermodynamic 
limit by using data for the largest two sizes (i.e., L — 15 and 
18). 



Ji jC is more prominent as 9 gets closer to 7r/4. This re- 
flects the fact that quantum fluctuations for the ID sys- 
tems become larger as 9 approaches to the SU(3) point, 
where more low-energy excitations appear. Except for 
the case of 9 = 0.3n, in which size effect may be pro- 
found, a smooth extrapolation of Ji. c to zero in the ther- 
modynamic limit is found for all 0's. This indicates that 
the 2D three-sublattice order will emerge for infinites- 
imal J i within the whole region of7r/4 < 9 < tt/2. 
In other words, the phase transitions out of the three- 
sublattice states actually occur at infinite a for original 
2D anisotropic BLBQ model. Thus our ED results lend 
strong support on the conclusions based on the scaling 
arguments discussed in the previous subsection. 
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IV. CONCLUSIONS 

To summarize, we elaborate the effect of spatial 
anisotropy a on the stability of the three-sublattice spin 
nematic state in the model of Eq. (1) through various 
analytic as well as numerical approaches. We conclude 
that the three-sublattice state is stable for all < a < oo 
within the whole region of ir/4 < 9 < ir/2. Our analysis 
thus gives a complete picture for smooth crossover from 
the triangular-lattice case to both the square-lattice and 
the ID-chain limits as the anisotropy a is varied. More- 
over, our work provides some insights on the validity of 
the LFW theory. Basically, the strength of the stability 
for the considered order can be understood within the 
simple LFW theory. Nevertheless, the predicted phase 
boundaries and the excitation spectrums is merely sug- 



gestive, especially in both of the a = and a->oo limits. 
Because the local constraints are released in the LFW 
analysis, it is interesting to see if great improvement can 
be obtained through other approaches (say, the varia- 
tional Monte Carlo method), in which these constraints 
are taken into account rigorously. Such discussions go be- 
yond the scope of the present work and deserve further 
investigations. 
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